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ABSTRACT 


This  report  presents  the  power-normal  and  -lognormal  models,  which  describe  the  effect  of 
specimen  size  on  the  distribution  of  life  or  strength  of  a product  or  material.  Such  a model  arises  when 
any  specimen  can  be  regarded  as  a series  system  of  smaller  portions,  where  portions  of  a certain  size  have 
a normal  or  lognormal  life  (or  strength)  distribution.  Also,  this  report  documents  the  first  computer 
program  that  fits  this  model  to  data  (including  censored  and  interval  life  data)  from  specimens  of  various 
sizes.  The  program  employs  maximum  likelihood  fitting  and  provides  approximate  confidence  limits,  as 
well  as  estimates,  for  model  parameters,  distribution  percentiles,  and  other  quantities  of  interest.  How 
to  run  the  program  is  explained  with  an  analysis  of  data  on  time  to  electromigration  failure  of  aluminum 
conductors  for  microcircuits.  The  program  runs  on  the  VAX  computer  and  SUN  workstation. 


1.  INTRODUCTION 

Purpose.  This  report  describes  in  detail  properties  of  the  power-lognormal  and  related  power- 
normal  models.  The  report  also  presents  the  first  computer  program  for  fitting  these  models  to  data.  This 
introduction  briefly  presents  the  history  of  the  models,  the  models,  the  computer  program,  and  an 
overview  of  this  report.  Then  the  introduction  describes  electromigration  failure  data  used  to  illustrate 
1)  use  of  the  models  of  Section  3 and  4 and  2)  how  to  run  the  computer  program  in  Section  2. 

1.1.  Background 

History.  The  power-normal  and  -lognormal  models  have  long  been  used  to  represent  the  effect 
of  specimen  size  on  the  distribution  of  life  or  strength  of  products  and  materials.  In  an  early  application, 
Chaplin  (1880)  used  the  power-normal  model  to  describe  the  distribution  of  tensile  strength  of  steel  bars 
of  various  lengths.  Nelson  (1990,  pp.  385-392)  presents  models  for  size  effect  and  applications  to  cable. 
Harter  (1977)  surveys  the  literature  on  this  and  other  models  for  the  effect  of  specimen  size  for  a great 
variety  of  applications.  Our  work  on  this  model  and  computer  program  was  motivated  by  an  application 
to  data  on  time  to  electromigration  failure  of  aluminum  conductors  in  microcircuits.  A key 
electromigration  reference  is  LaCombe  and  Parks  (1986),  who  use  the  model  for  their  data. 

The  models.  The  power-normal  and  power-lognormal  models  are  briefly  described  here.  Full 
details  appear  in  Sections  3 and  4.  For  specimens  of  size  S , the  reliability  function  of  the  power-normal 
distribution  is  the  population  fraction  above  the  value  y (of  life,  strength,  etc.);  namely, 

F(y,S;n,a,S J = {$t-(y-^)/a]}5/5°,  -oo  < y < ».  (1.1) 

This  is  a power  of  the  normal  reliability  function.  Here  <i>[  ] is  the  standard  normal  cumulative 
distribution  function  and  the  minus  sign  converts  it  into  a reliability  function.  The  location  parameter  \i 
has  the  range  — oo  < ^ < » . The  scale  parameter  a is  positive.  The  normal  size  parameter  S0  is  positive; 
that  is,  y has  a normal  distribution  when  S = S0.  In  previous  work,  people  had  to  assume  a value  for  So-> 
which  is  unknown.  The  computer  program  here  provides  an  estimate  of  S0  (and  of  p.  and  a)  from  data. 
The  ratio  p = S/S0  is  called  the  power  parameter.  The  reliability  function  of  the  power-lognormal 
distribution  is  the  population  fraction  above  the  value  t;  namely, 
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(1.2) 


Git,S;n,<y,SJ  = {^-(ln(r)-M)/a]}5/5°,  t > 0. 

The  parameters  n and  a have  the  same  names  and  ranges.  Here  S0  is  the  lognormal  size  parameter ; that 
is,  t has  a lognormal  distribution  when  S = S0. 

Computer  program.  The  computer  program  POWNOR  uses  maximum  likelihood  (ML)  methods 
to  fit  the  model  to  data.  Program  output  includes  ML  estimates  and  approximate  confidence  limits  for 
model  parameters,  distribution  percentiles,  and  other  quantities  of  interest.  User  friendly,  the  program 
is  written  in  Fortran  77  and  runs  on  the  VAX  1 1/785  computer  under  the  VMS  operating  system  and  on 
the  SUN  workstation.  It  employs  IMSL  routines.  Copies  of  the  program  may  be  requested  from  the 
Statistical  Engineering  Division,  Adm.  Bldg.  A337,  National  Institute  of  Standards  and  Technology, 
Gaithersburg,  MD  20899,  (301)  975-2839. 

Overview.  This  report  contains: 

• Section  2 on  how  to  run  the  computer  program,  illustrated  with  data  on  time  to  electromigration 

failure.  Many  readers  need  to  read  only  Section  2. 

• Section  3 on  the  statistical  properties  of  the  power-normal  model. 

• Section  4 on  the  statistical  properties  of  the  related  power-lognormal  model. 

• Section  5 on  the  maximum  likelihood  theory  for  fitting  either  model  to  data,  which  may  be 

censored  (unfailed  runouts)  or  interval  data. 

• Section  6 on  the  numerical  methods  employed  to  carry  out  the  maximum  likelihood  calculations. 

Acknowledgments.  The  authors  gratefully  acknowledge  the  support  of  this  research  under  their 
NIST/NSF/ASA  Research  Fellowship  at  the  National  Institute  of  Standards  and  Technology  (NIST).  The 
National  Science  Foundation  (NSF)  generously  provided  our  salary  and  travel  expenses.  NIST  kindly 
provided  research  facilities  and  colleagues.  The  American  Statistical  Association  (ASA)  handled  the 
Fellowship  selection  and  administration.  Dr.  Robert  Lundegard  and  Mrs.  Ruth  Varner  of  NIST 
effectively  and  graciously  administered  our  work  at  NIST.  Mr.  Harry  Schafft  and  Dr.  Jim  Lechner  of 
NIST  provided  much  valued  stimulation  for  this  report  in  our  collaboration  on  electromigration  research 
with  them. 

1.2  Example  Data 

Overview.  This  section  presents  electromigration  failure  data  used  to  illustrate  how  to  run  the 
program  to  fit  a power-(log)normal  model  to  data.  The  section  first  gives  background  on 
electromigration.  Then  it  describes  the  data.  Finally  it  states  a purpose  of  fitting  a power-lognormal 
model  to  the  data.  The  results  of  fitting  appear  in  Section  2. 

Electromigration.  Aluminum  conductors  in  microcircuits  can  fail  from  movement  of  atoms 
resulting  in  an  open  conductor  or  in  a short  (due  to  an  aluminum  filament  growing  between  adjoining 
conductors).  Such  failures  occur  more  quickly  at  high  current  density  and  high  temperature.  Microcircuit 
components  and  conductors  are  being  built  smaller  and  smaller  for  higher  speed,  resulting  in  higher 
current  densities  and  temperatures  in  conductors.  The  microcircuit  industry  fears  that  electromigration 
failure  may  limit  such  size  reduction.  So  the  industry  seeks  more  knowledge  of  electromigration  including 
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a model  for  the  distribution  of  time  to  failure  and  for  the  effect  of  conductor  length.  The  power- 
lognormal  model  has  been  proposed,  but  some  experts  have  observed  that  such  time  to  failure  data  for 
various  specimen  lengths  all  appear  adequately  described  by  a lognormal  distribution.  The  following  data 
yield  information  on  the  adequacy  of  the  power-lognormal  and  lognormal  distributions. 

Data.  Table  1.1  displays  hours  to  failure  of  59  test  conductors  of  length  S = 400  microns.  All 
specimens  ran  to  failure  at  a certain  high  temperature  and  current  density.  The  data  come  from  an 
interlaboratory  comparison  reported  by  Schafft  and  others  (1987).  The  59  specimens  were  all  tested  under 
the  same  temperature  and  current  density.  The  data  appear  in  a lognormal  probability  plot  in  Figure  1.1. 
The  plot  is  relatively  straight.  Thus  it  suggests  that  a lognormal  distribution  adequately  fits  the  data  and 
that  a power-lognormal  distribution  fitted  to  the  data  is  "close"  to  lognormal;  that  is,  p = S/S0  is  "close" 
to  1. 


Table  1.1  Hours  to  Failure  of  400-micron  Electromigration  Specimens 


6.545 

9.289 

7.543 

5.956 

6.492 

5 .459 

S . 120 

4 .706 

8.637 

2 . 997 

8.591 

6.129 

11. 033 

5.331 

6 .953 

4 . 283 

6 . 522 

4.137 

7 .459 

7.495 

6.573 

6 . 533 

5 . 539 

6.037 

5.307 

5.725 

8.532 

9.563 

6.369 

7.024 

3.336 

9 .218 

7.945 

6 .869 

6.352 

4 .700 

6 . 54S 

9 . 254 

5 . 009 

7.439 

7 .393 

6.033 

10 . 092 

7.496 

4.531 

7 . 974 

3 .799 

7.633 

7.224 

6.515 

7.365 

6.476 

6 . 923 

6 . 071 

5 . 540 

10 . 491 

5.434 

5 .923 

7.937 

tlM 


Fig.  1.1  Lognormal  Probability  plot  of  400-micron  data. 
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Purpose.  Fitting  the  power-lognormal  distribution  to  the  data  serves  certain  purposes.  First,  it 
provides  an  estimate  of  the  lognormal  length  which  previously  has  not  been  properly  estimated  and 
may  lead  to  physical  insight.  Second,  wide  confidence  limits  for  p = S/S0  that  enclose  1 (the  lognormal 
distribution)  would  indicate  why  many  data  sets  appear  adequately  described  by  a lognormal  distribution. 
Third,  extrapolation  to  estimate  the  life  distribution  of  the  total  length  of  conductor  in  microcircuits  under 
normal  temperature  and  current  density  is  of  interest.  Fourth,  the  uncertainty  of  that  distribution  estimate 
provides  guidance  on  the  appropriate  number  of  test  specimens  needed  to  yield  estimates  of  a desired 
accuracy. 


2.  HOW  TO  RUN  THE  POWNOR  PROGRAM 

Introduction.  This  section  explains  how  to  run  the  program  POWNOR  to  fit  a power-normal 
or  -lognormal  model  to  data.  It  describes  and  illustrates  input  files,  program  operation,  and  printed 
output.  The  examples  here  use  the  electromigration  data  of  Section  1.2. 

Input  files.  The  user  must  create  two  input  files  before  the  program  can  be  run.  One  of  these 
files  POWNOR. DATA  stores  data  on  specimen  failure  times  and  specimen  lengths.  The  other  file 
POWNOR. OPTIONS  contains  various  user  choices  (e.g.,  desired  level  of  detail  in  printing  results, 
starting  values  for  parameters,  etc.).  The  user  can  employ  his  favorite  text  editor  to  create  the  input  files. 
In  both  of  these  files,  values  in  a line  are  separated  by  one  or  more  spaces  (free  format).  Each  line 
contains  at  most  80  characters  and  spaces  (typical  line  length  on  many  computer  terminals). 

Data  file.  This  file  contains  data  on  specimen  failure  times  and  specimen  lengths.  Figure  2.1 
shows  the  electromigration  data  as  they  appear  in  the  data  file  POWNOR.DATA.  On  each  line  (record) 
of  this  file,  enter  the  data  on  a particular  specimen.  The  three  values  on  the  line  for  a specimen  are: 

1 . The  lower  endpoint  of  the  censoring  interval.  This  data  value  is  numerical  and  must  not  be  below 
— IE  10. 

2.  The  upper  endpoint  of  the  censoring  interval.  This  data  value  is  numerical  and  must  not  exceed 
1E10. 

3.  Length  of  the  specimen.  This  data  value  is  numerical  and  cannot  be  zero  or  negative. 


Censoring  is  indicated  as  follows: 

• If  the  specimen  fails  at  an  exactly  known  age,  set  the  two  endpoints  of  the  censoring  interval 
equal  to  the  failure  age. 

• If  the  failure  is  in  a known  interval,  set  the  lower  and  upper  endpoints  equal  to  the  interval 
endpoints. 

• If  the  specimen  is  left  censored,  set  the  lower  endpoint  equal  to  — 1E10,  and  the  upper  value 
equal  to  the  censoring  age. 

• If  the  specimen  is  right  censored,  set  the  lower  endpoint  equal  to  the  censoring  age  and  the  upper 
endpoint  equal  to  IE  10. 
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LOWER 

6 . 545 
8.120 
11.038 

7.459 
5.807 
8 .336 
6.948 

10.092 

7.224 

6.515 

9.289 

4.706 

5.381 

7.495 
6.725 
9.218 
9.254 

7.496 
7.365 
6 . 476 
7.543 
8 . 637 
6.958 
6.573 
8 . 532 
7.945 
5.009 
4.531 

6.923 
6.071 
6.956 
2.997 
4.288 
6.538 
9.663 
6.869 
7.489 
7.974 
5.640 

10.491 
6.492 
8.591 
6.522 
5.589 
6.369 
6.352 
7.398 
8.799 
5 .434 

5.923 

5.459 
6 . 129 
4 . 137 
6 .087 
7 . 024 
4 .700 
6.033 
7 . 683 
7.937 


UPPER 

LENGTH 

6.545 

400 

8 . 120 

400 

11.038 

400 

7 .459 

400 

5.807 

400 

8 .336 

400 

6.948 

400 

10 . 092 

400 

7 .224 

400 

6.515 

400 

9 .289 

400 

4.706 

400 

5.381 

400 

7.495 

400 

6.725 

400 

9.218 

400 

9.254 

400 

7 . 496 

400 

7.365 

400 

6 . 476 

400 

7.543 

400 

8 . 687 

400 

6.958 

400 

6.573 

400 

8.532 

400 

7.945 

400 

5 . 009 

400 

4.531 

400 

6 . 923 

400 

6.071 

400 

6.956 

400 

2.997 

400 

4.288 

400 

6.538 

400 

9.663 

400 

6.869 

400 

7.489 

400 

7.974 

400 

5.640 

400 

10.491 

400 

6.492 

400 

8.591 

400 

6.522 

400 

5.589 

400 

6 .369 

400 

6.352 

400 

7.398 

400 

8.799 

400 

5.434 

400 

5.923 

400 

5 .459 

400 

6.129 

400 

4 . 137 

400 

6.087 

400 

7.024 

400 

4 .700 

400 

6.033 

400 

7 . 633 

400 

7 .937 

400 

Fig.  2.1  Electromigration  life  data  in  the  file  POWNOR.DATA. 
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Section  5.2  defines  different  types  of  censored  data.  The  electromigration  data  in  Figure  2.1  do  not 
include  any  censored  observations.  In  this  example,  each  specimen  has  the  same  length,  400  microns. 
The  file  must  not  contain  the  headings  LOWER,  UPPER,  and  LENGTH  added  to  Figure  2.1  for 
explanation.  The  data  file  cannot  contain  more  than  10,000  lines. 

Options  file.  The  POWNOR. OPTIONS  file  contains  the  following  choices  of  the  user,  which 
are  explained  later: 

• Whether  or  not  log  transformation  of  failure  data  in  POWNOR. DATA  is  required. 

• Level  of  detail  in  documenting  the  iteration  history. 

• Starting  values  for  parameters  to  be  estimated  and  fixed  parameter  values. 

• Specimen  lengths  for  which  estimates  of  distribution  percentiles  and  their  associated  confidence 
limits  are  desired. 

Figure  2.2  shows  the  POWNOR. OPTIONS  file  used  for  the  electromigration  data  example. 


1 

1 

0 2.9219 

0 -0.7186 

0 2.5415 

400 


Fig.  2.2  The  POWNOR. OPTIONS  file  with  options  for  the  electromigration  data. 


Log  transformation.  Enter  1 on  the  first  line  to  calculate  natural  logs  of  the  failure  (and 
censoring)  times  in  POWNOR. DATA.  If  no  log  transformation  is  required,  enter  0.  The 
electromigration  example  requires  log  transformation  of  failure  times.  This  is  provided  by  the  value  1 
on  the  first  line  of  the  file  POWNOR. OPTIONS  in  Figure  2.2.  The  program  fits  the  power -normal 
distribution  to  the  log  data. 

Iteration  history.  Enter  1 on  the  second  line  for  detailed  printout  of  iterations.  The  program 
then  stores  the  iterations  output  in  a file  named  POWNOR. ITERATIONS.  Enter  0,  if  the  detailed 
iteration  history  is  not  required.  This  output  is  recommended  as  an  aid  to  understanding  iterative  fitting 
that  is  not  converging  well.  The  POWNOR. OPTIONS  file  in  Figure  2.2  requests  that  the  detailed 
iterations  be  saved. 

Parameter  values.  The  user  is  required  to  supply  starting  values  of  the  parameters  for  numerical 
maximization  of  the  sample  log  likelihood  function.  Section  6 provides  guidelines  to  obtain  starting 
values.  Alternatively,  the  user  may  wish  to  hold  some  of  the  parameters  fixed  at  desired  values.  The 
lines  3,  4,  and  5 of  POWNOR. OPTIONS  specify  such  choices  for  /z,  ln(cr),  and  ln(S0),  respectively.  Each 
of  these  lines  must  contain  2 values.  A 0 for  the  first  value  indicates  that  the  parameter  is  to  be  estimated 
from  data.  The  second  value  is  then  taken  to  be  the  starting  value  for  the  parameter.  A 1 for  the  first 
value  indicates  that  the  parameter  value  is  held  fixed  at  the  second  value.  This  feature  is  particularly 
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useful  in  likelihood  ratio  tests  for  parameters,  as  described  in  Nelson  (1990,  Chapter  9).  Figure  2.2 
indicates  that  all  three  parameters  in  the  electromigration  data  example  are  to  be  estimated  from  data. 
The  starting  values  for  the  /x,  ln(a),  and  ln(50)  are  2.9219,  —0.7186,  and  2.5415,  respectively.  These 
starting  values  were  obtained  following  the  guidelines  in  Section  6. 

Specimen  lengths  for  percentile  estimates.  Line  6 and  the  ones  below  it  contain  the  specimen 
lengths  for  which  estimates  and  confidence  limits  for  distribution  percentiles  are  desired.  Each  line 
contains  a single  length  value.  The  user  can  specify  at  most  10  lengths.  In  the  electromigration  example, 
such  percentile  estimates  are  requested  only  for  specimens  of  length  400  microns. 

Program  operation.  The  user  must  previously  create  the  two  input  files  as  described  above.  The 
program  then  can  be  run  by  typing  RUN  POWNOR  at  the  $ prompt  of  a VAX  computer  with  a VMS 
operating  system.  Figure  2.3  shows  program  messages  during  run  for  the  electromigration  data  example. 
The  underlined  entries  in  the  figure  are  user  responses. 


$ RUN  POWNOR 

PROGRAM  MESSAGE:  The  program  POWNOR  reads  the 

user  options  from  POWNOR.  OPTIONS 
and  data  lines  from  POWNOR. DATA 

PROGRAM  MESSAGE:  Iteration  history  is  in  POWNOR. ITERATIONS 

PROGRAM  MESSAGE:  Program  output  is  in  POWNOR. OUTPUT 

FORTRAN  STOP 

$ TYPE  POWNOR  .OUTPUT 


Fig.  2.3  How  to  submit  a POWNOR  run  on  a VAX  computer. 


Output.  Figure  2.4  displays  the  program  output  in  the  file  POWNOR. OUTPUT  from  the  ML 
fitting  of  a power-normal  model  to  the  logs  of  the  electromigration  data.  This  output  is  explained  below; 
the  numbers  correspond  to  those  in  the  figure. 

I.  This  output  states  the  number  of  uncensored  and  censored  observations  in  the  data  file 
POWNOR.  DATA. 

II.  This  output  shows  the  ML  estimates  of  the  (transformed)  model  parameters  Ox,lncr,lnS0)  and  their 
approximate  95%  confidence  limits.  For  example,  the  estimate  of  a is  o = exp(— 0.7336)  = 
0.48,  and  its  95%  confidence  limits  are  a = exp(-1.6375)  = 0.19  and  a = exp(0.1704)  = 
1.19. 
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POWNOR  Version  1.0 


3 O-AUG-9 1 11:50:49 


Developed  by  Necip  Doganaksoy  and  Wayne  Nelson  under  1991  Fellowship 
grant  from  ASA/NSF/NIST.  The  program  is  documented  in  NIST  Report#. 


Il  Total  number  of  data  cases=  59 

Exact  observations  = 59 

Left  censored  observations=  0 


Right  censored  observations  = 0 
Interval  censored  observations=  0 


Maximized  log-likelihood  = 55.8554 

Maximum  likelihood  estimates  for  distribution  parameters 
with  approximate  95%  confidence  limits 


II: 


III: 


IV: 


Parameter 

ML  Estimate 

Lower  Limit 

Upper  Limit 

1 mu 

2.8779 

0.8280 

4 . 9279 

2 ln(sigma) 

-0.7336 

-1.6375 

0.1704 

3 In (normal  length) 

2.6779 

-3 . 1808 

8 . 5366 

Estimated 

information  matrix 

1 

2 

3 

1 

1159.60 

-985.06 

255.65 

2 

-985.06 

947.61 

-200.44 

3 

255.65 

-200.44 

59.00 

Estimated 

covariance  matrix 

1 

2 

3 

1 

1.09391 

0.47809 

3.11576 

2 

0.47809 

0.21270 

1.34900 

3 

-3 . 11576 

-1.34900 

8.93479 

Estimated 

correlation  matrix 

1 

2 

3 

1 

1.00000 

0.99115 

0.99662 

2 

0.99115 

1.00000 

0.97856 

3 

-0.99662 

-0.97856 

1.00000 

Maximum 

likelihood  estimates  for  distribution  percentiles 

with  approximate  95% 

confidence  limits 

Length= 

400.0000 

Pet. 

ML  Estimate 

Lower  Limit 

Upper  Limit 

Std.  Error 

0.1 

0.9731 

0.5994 

1.3463 

0.1907 

0.2 

1.0541 

0.7256 

1.3825 

0.1676 

0.5 

1.1664 

0.8968 

1.4360 

0.1375 

1 

1.2561 

1.0297 

1.4825 

0.1155 

2 

1.3507 

1.1655 

1.5358 

0.0945 

5 

1.4853 

1.3496 

1.6210 

0.0692 

10 

1.5973 

1.4930 

1.7016 

0.0532 

30 

1.8074 

1.7357 

1.8792 

0.0366 

50 

1.9360 

1.8719 

2 . 0001 

0.0327 

70 

2 . 0528 

1.9921 

2 . 1135 

0.0310 

90 

2 . 2044 

2 . 1380 

2 . 2709 

0.0339 

95 

2.2714 

2 . 1944 

2 . 3485 

0 . 0393 

98 

2 .3431 

2.2475 

2 .4387 

0.0488 

99 

2 .3890 

2.2776 

2 . 5004 

0.0568 

99 . 5 

2 .4298 

2 .3020 

2 . 5576 

0.0652 

99 . 9 

2 .5108 

2 . 3439 

2 . 6776 

0.0851 

Std. 


Fig.  2.4  POWNOR  output  for  the  electromigration  data. 


Error 
1.0459 
0.4612 
2 . 9891 
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III.  The  local  estimates  of  the  information,  covariance,  and  correlation  matrices  for  the  (transformed) 
parameter  estimates  are  shown  here.  The  high  correlations  (near  + 1)  indicate  that  convergence 
to  the  ML  estimates  may  be  incomplete  and  that  a better  parameterization  of  the  model  could 
avoid  convergence  problems. 

IV.  This  output  shows  the  ML  estimates  of  percentiles  (log  values)  and  their  95%  confidence  limits 
at  each  specified  length. 

The  theory  underlying  the  ML  fitting  is  in  Section  5.  Section  6 explains  numerical  methods  used  in 
fitting. 


Reparametrization.  The  program  POWNOR  maximizes  the  log  likelihood  function  with  respect 
to  p,  ln(a),  and  ln(S0).  Sometimes,  this  parametrization  of  the  model  must  be  altered  to  improve 
numerical  convergence  to  the  ML  estimates.  Linear  orthogonal ization  of  the  parameters  provides  a 
convenient  means  for  improving  convergence  as  described,  for  example,  in  Ross  (1990).  Section  6 
describes  such  linear  orthogonalization  of  parameters.  Another  version  of  the  program 
POWNOR_ORTHO  uses  such  orthogonalization.  The  input  files,  operation,  and  outputs  of 
POWNOR_ORTHO  are  very  similar  to  those  of  POWNOR. 


3.  THE  POWER-NORMAL  DISTRIBUTION 


Introduction.  This  section  presents  properties  of  the  power-normal  distribution.  As  background, 
readers  need  to  know  basic  properties  of  the  normal  distribution  and  normal  probability  paper.  For 
example.  Nelson  (1982,  1990)  covers  these  topics  and  the  series-system  model  (below)  in  detail.  The 
following  topics  include  its  reliability  function,  cumulative  distribution  function,  the  series-system  model, 
percentiles,  median,  probability  density,  mode,  mean,  variance,  standard  deviation,  and  hazard  function 
(failure  rate). 


Reliability  function  and  cdf.  The  reliability  Junction  of  the  power-normal  distribution  is  the 
population  fraction  above  the  value  y;  namely. 


F(y;p,a,p)  ■ {${-(y-p)/o]}‘) , -oo  < y < oo. 


(3.1) 


Here  - oo  < p < oo  is  the  location  parameter,  a > 0 is  the  scale  parameter,  and  p > 0 is  the  power 
parameter  p = S/S0  in  (1.1).  Note  the  minus  sign  in  the  standard  normal  cumulative  distribution  function 
$[  ];  it  yields  the  upper  tail  probability  of  the  normal  distribution.  For  p = 1,  the  distribution  is  normal; 
only  then  is  p the  mean  and  median  of  the  distribution,  and  only  then  is  a the  standard  deviation  of  the 
distribution.  The  cumulative  distribution  function  (cdf)  is  the  population  fraction  below  y;  namely, 

F(y;p,o,p ) = 1 - {<t(-(y-/i)/a]}p,  -oo  < y < oo.  (3-2) 

The  distribution  is  simpler  in  terms  of  the  standardized  deviate  z = (y-p)/a.  Then 

F(y;p,o,p ) = [<$(-z)]p,  F(y;p,o,p)\  = 1 - [<i<-z)]p,  -oo  < z < oo.  (3-3) 


David  (1981,  p287)  lists  tables  of  these  functions  for  integer  p values. 
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Fig.  3.1  Power-normal  cdfs  on  normal  probability  paper. 


Probability  plot.  Figure  3.1  displays  power-normal  cdfs  for  selected  p values  on  normal 
probability  paper.  For  p = 1,  the  cdf  is  normal  and  is  a straight  line.  For  p > 1,  the  cdfs  curve 
concave  to  the  left  and  shift  toward  lower  y (and  z ) values  as  p gets  larger.  For  0 < p < 1,  the  cdfs 
curve  concave  to  the  right  and  shift  toward  higher  y (and  z ) values  as  p gets  smaller. 

Approximate  cdf.  For  large  p,  the  exact  cdf  formula  (3.2)  requires  a table  or  algorithm  for  <£(  ) 
that  is  accurate  in  the  far  tails.  Asymptotically  as  p -*  oo,  the  power-normal  cdf  approaches  the  cdf  of 
a smallest  extreme  value  distribution  (convergence  in  law  or  in  distribution).  Bury  (1975)  gives  the 
location  and  scale  parameters  of  that  distribution  as  a function  of  p,  a,  and  p.  Convergence  to  this 
asymptotic  distribution  is  slow.  That  is,  p must  be  very  large  (say  well  above  100)  for  a satisfactory 
approximation. 

Series-system  model.  In  some  applications,  product  comes  in  different  sizes,  and  its  life, 
strength,  or  some  property  depends  on  size.  For  example,  time  to  electromigration  failure  of  conductors 
in  microcircuitry  depends  on  conductor  length.  Also,  strength  of  metal  bars  depends  on  bar  length 
(Chaplin  (1880)).  Harter  (1977)  surveys  models  for  the  effect  of  size.  The  following  series-system  model 
is  a simple,  widely-used  model  for  the  effect  of  size.  For  purposes  of  presenting  the  model,  reliability 
terminology  and  product  life  are  used.  The  assumptions  of  the  series-system  model  for  size  effect  are: 

1 . If  product  is  conceptually  divided  into  portions  of  any  size  of  interest,  it  is  assumed  that  the  times 
to  failure  of  the  portions  are  statistically  independent. 

2.  The  product  fails  when  the  first  such  portion  fails. 

3.  Portions  of  size  S0  have  a reliability  function  R0(t),  which  is  the  population  fraction  surviving 
beyond  age  t. 
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(3.4) 


It  follows  from  these  assumptions  that  product  of  size  S has  a reliability  function 

Rs(t)  = [ff0(/)]“0. 

Nelson  (1982,  1990)  presents  this  model  in  more  detail  and  gives  various  applications.  When  R0(t ) is  the 
reliability  function  of  a normal  distribution,  (3.4)  is  the  power-normal  reliability  function  (3.1).  Then  the 
power  parameter  p = S/S0  is  the  ratio  of  the  product  size  S to  the  "normal  size"  S0. 


Fractiles.  The  F tractile  of  the  power-normal  distribution  is 

VF  = n * zF-a;  (3.5) 

here  zr  is  the  F'  = 1 - (1  -F)I/p  fractile  of  the  standard  normal  distribution.  For  p > > 1,  F'  usually 
is  very  near  zero.  For  p near  zero,  F'  usually  is  very  near  one.  Tables  of  zF>  appear  in  most  statistics 
books.  However,  few  tables  adequately  provide  zF>  values  for  F'  very  near  zero  or  one. 

Median.  The  median  is  the  0.50  fractile  7 j0  50.  It  is  a "center"  of  the  distribution  in  that  half  of 
the  population  is  below  77  50  and  half  is  above  it.  For  the  power-normal  distribution, 

^0.50  = M + zr (3>6) 

here  F’  = 1 — (0.5)1/p.  For  p = 1,  the  distribution  is  normal,  F = 0.50,  zF  = 0,  and  tj0  50  = p. 


Density.  The  power-normal  probability  density  is  the  derivative  of  the  cdf  (3.2);  namely, 

fiy;n,t.p)  - (pic)  m -p)lc]  W-ly-tilc}]"-'-,  (3  7) 

here  0(z)  = (27r)_1/2  exp(-z2/2)  is  the  standard  normal  probability  density.  In  terms  of  the  standardized 
deviate  z = (y-p)/a, 

fiy;p,a,p)  = (pic ) tfz)  [^-z)]"'1.  (3-8> 

Figure  3.2  displays  such  probability  densities  as  a function  of  z for  selected  p values  where  p = 0 and 
cr  = 1.  For  p = 1,  the  density  is  the  standard  normal  density.  These  densities  are  unimodal. 

SHO-tOOQO  3000  1000  300  100  *0  20  3 1 0.5  0.2  0.1 


2 

Fig.  3.2  Power-normal  probability  densities. 
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Mode.  The  mode  of  a power-normal  distribution  is  the  value  i]m  where  the  probability  density 
peaks.  It  is  the  most  probable  y value.  t]m  is  the  solution  of  dfly;p,o,p)/dy  = 0;  namely, 

1m  ' /*  + Zm0'  <39) 

Here  zm  is  a function  of  just  p and  is  the  solution  of  -zm$(-zm)/<t>(zm)  - p — 1.  For  p = 1,  zm  = 0 and 
the  mode  r;m  = p,  a well-known  property  of  the  normal  distribution. 

Mean.  The  mean  of  the  power-normal  distribution  is 

EY  = \ yfy;p,a,p ) dy  = p + Epa.  (3-10) 

Here  the  mean  Ep  for  the  standardized  distribution  (p  = 0 and  a — 1)  is  a function  of  p only;  namely, 

£p  ■ l T.w  #z>  w-z)]''1  *•  (3-U) 

For  integer  p,  Ep  is  the  mean  of  the  first  order  statistic  of  a sample  of  size  p.  David  (1981,  p291)  lists 
tables  of  Ep  for  integer  p.  Ep  apparently  has  not  been  tabulated  for  non-integer  p;  then  (3.11)  must  be 
evaluated  by  numerical  quadrature. 

Variance  and  standard  deviation.  The  variance  of  the  power-normal  distribution  is 

V(Y)  - | °‘jy-EY?fly,w,p)dy  = V /■  (3.12) 

here  the  variance  Vp  for  the  standardized  distribution  (p  = 0 and  a = 1)  is 

v,  = |_"<>(z-£p2<><«z)W-z)r1&.  (3-13) 

For  integer  p,  Vp  is  the  variance  of  the  first  order  statistic  of  a sample  of  size  p from  a standard  normal 

distribution.  David  (1981,  p291)  lists  tables  of  Vp  for  integer  p.  Vp  apparently  has  not  been  tabulated 

for  non-integer  p;  then  (3.13)  must  be  evaluated  by  numerical  quadrature.  The  standard  deviation  of  the 
power-normal  distribution  is 

sd{Y)  = [V{Y)]m  = (3-14) 

For  p = 1,  the  distribution  is  normal,  and  Vj  = 1,  V(Y)  — a2,  and  sd(Y)  = a. 

Hazard  function.  Much  used  in  reliability  work,  the  hazard  function  (also  called  the 
instantaneous  failure  rate ) of  the  power-normal  distribution  is 

h(y;p,o,p)  m fy;p,a,p)/[l-F(y;p,o,p)] 

= (p/o)<t>[(y-p)/o]/${-(y-p)/o] 

= p-h[(y-n)la]  = ph(z).  (3'15) 

Here  h(z)  = 4>(z)/$(—z)  is  the  hazard  function  of  the  standard  normal  distribution  and  is  an  increasing 
function  of  z = (y—p)lo.  Figure  3.3  depicts  h(z).  Any  power-normal  hazard  function  is  merely  a 
multiple  p of  h{z). 
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Fig.  3.3  Standard  power-normal  hazard  function. 


4.  THE  POWER-LOGNORMAL  DISTRIBUTION 


Introduction.  This  section  presents  properties  of  the  power-lognormal  distribution.  As 
background,  readers  need  to  know  basic  properties  of  the  power-normal  and  normal  distributions  and 
lognormal  probability  paper.  For  example,  Nelson  (1982,  1990)  covers  these  topics  and  the  series-system 
model  (below)  in  detail.  The  following  topics  include  its  reliability  and  cumulative  distribution  functions, 
its  relationship  to  the  power-normal  distribution,  the  series-system  model,  percentiles,  median,  probability 
density,  mode,  mean,  variance,  standard  deviation,  and  hazard  function  (failure  rate). 

Reliability  function.  The  reliability  function  of  the  power-lognormal  distribution  is  the  population 
fraction  surviving  beyond  age  t > 0;  namely, 

- {4(-(lii(f)-(*)/<7]}'.f>0.  (4-[) 

Here  <i>[  ] is  the  standard  normal  cumulative  distribution  function,  and  p > 0 is  called  the  power 
parameter.  Also  p (-<»  < p < co)is  called  the  location  parameter,  but  is  better  called  "mu".  Strictly 
speaking,  p is  a location  parameter  only  for  the  power -normal  distribution  when  y = ln(f)  is  used.  For 
the  power-/pgnormal  distribution,  r = e^  is  a scale  parameter  as  shown  below.  Similarly,  a > 0 is 
called  the  scale  parameter,  but  is  better  called  "sigma".  Strictly  speaking,  a is  a scale  parameter  only 
for  the  power -normal  distribution,  when  y = ln(r)  is  used.  For  the  power-/ognormal  distribution,  p and 
a are  actually  shape  parameters.  For  p = 1,  the  distribution  is  lognormal;  only  then  is  p the  mean  of  log 
life  y = ln(r),  is  the  median  life,  and  is  a the  standard  deviation  of  log  life  y = ln(r).  The  lognormal 
distribution  is  briefly  presented  in  Nelson  (1982,  1990)  and  in  many  other  statistics  books;  Shimizu  and 
Crow  (1988)  devote  a book  to  it. 


Cdf. 

namely, 


The  cumulative  distribution  function  (cdf)  is  the  population  fraction  failing  before  age  r; 


G(r,p,o,p ) = 1 - (f{-(ln(r)-p)/a]}p,  t > 0. 


(4.2) 
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The  distribution  is  simpler  in  terms  of  the  log-transformed  standardized  deviate  z = [ln(f)— /z]/ct.  Then 

G(r,p,o,p)  - [4<-z)]»,  G(t;p,o,p)  - 1 - [<K-z)]<\  (4'3) 

David  (1981,  p287)  list  tables  of  these  functions  for  integer  p values. 

Probability  plot.  Figure  4.1  displays  power-lognormal  cdfs  for  selected  p values  on  lognormal 
probability  paper.  Figure  4.1  is  the  same  as  Figure  3.1,  except  that  Figure  4.1  has  a data  scale  (for  t ) 
that  is  logarithmic  and  Figure  3. 1 has  a data  scale  (for  y)  that  is  linear.  For  p = 1,  the  cdf  is  lognormal 
and  is  a straight  line.  For  p > 1,  the  cdfs  curve  concave  to  the  left  and  shift  toward  lower  t (and  z ) 
values  as  p gets  larger.  For  0 < p < 1,  the  cdfs  are  curves  concave  to  the  right  and  shift  toward  higher 
t (and  z ) values  as  p gets  smaller.  Figure  4.2  shows  selected  cdfs  that  have  common  values  for  the  0.25 
and  0.75  fractiles  to  give  these  distributions  roughly  the  same  spread,  as  if  fitted  to  the  same  data  set. 
The  figure  shows  that  the  lower  tails  differ  appreciably.  Thus  it  is  important  to  know  p if  extrapolating 
into  the  lower  tail  in  high  reliability  applications. 

Other  parameters.  For  some  purposes,  it  is  useful  to  express  the  distribution  in  terms  of  its  true 
scale  parameter  r = exp (p).  Then 


C{t-p,a,p)  = {4f-(ln(//7))/o])p  = {4t-(In(f/r)'/»)J}p, 

G(t;p,a,p)  = 1 - {$(-(ln(//T))/a]}p  = 1 - {<S{-(ln(r/r)1'«)]}‘’. 
For  p = 1,  t is  the  distribution  median;  otherwise,  it  is  the  1 — (0.5)°  ffactile. 


(4.1’) 

(4.2') 


Approximate  cdf.  For  large  p,  the  exact  cdf  formula  (4.2)  requires  a table  or  algorithm  for  <£(  ) 
that  is  accurate  in  the  far  tails.  Asymptotically  as  p -*  oo,  the  power-lognormal  cdf  approaches  the  cdf 
of  a Weibull  distribution  (convergence  in  law  or  in  distribution).  Bury  (1975)  gives  the  scale  and  shape 
parameters  of  that  Weibull  distribution  as  a function  of  p,  a,  and  p.  Convergence  to  this  asymptotic 
Weibull  distribution  is  slow.  That  is,  p must  be  very  large  (say,  well  above  100)  for  a satisfactory 
approximation. 

Series-system  model.  As  described  in  Section  3,  a product  may  come  in  different  sizes,  and  its 
life  may  depend  on  size.  The  previous  assumptions  for  a series-system  model  are  assumed  here. 
However,  if  in  (3.4)  the  reliability  function  R0(t ) for  product  of  size  S0  is  lognormal,  then  the  reliability 
function  (3.4)  of  product  of  size  S is  the  power-lognormal  one  (4. 1).  Then  the  power  parameter  p = S/S0 
is  the  ratio  of  the  product  size  S to  the  "lognormal  size"  S0. 

Relation  to  the  power-normal.  The  power-lognormal  distribution  has  a relationship  to  the 
simpler  power-normal  distribution.  If  t has  a power-lognormal  distribution  with  parameters  p,  o,  and  p, 
then  y = ln(r)  has  a power-normal  distribution  with  the  same  values  of  the  parameters  p,  a,  and  p.  Thus, 
one  can  analyze  power-lognormal  data  by  fitting  the  power-normal  distribution  to  the  log  data  y = ln(r). 
The  power-normal  distribution  is  simpler  in  that  only  p is  a shape  parameter,  whereas,  for  the  power- 
lognormal  distribution,  both  p and  o are  shape  parameters. 


14 


PERCENT  FAILED  PERCENT  FAILED 


Fig.  4.1  Power-lognormal  cdfs  on  lognormal  probability  paper. 
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Fractiles.  The  F fractile  of  the  power-lognormal  distribution  is 

tF  = exp(/ i+zF-o);  (4.4) 

here  zr  is  the  F'  = 1 - (1  -F)llp  fractile  of  the  standard  normal  distribution.  tF  is  the  antilog  of  the  F 
fractile  (3.5)  of  the  corresponding  power-normal  distribution. 

Median,.  The  median  of  the  power-lognormal  distribution  is 

*o.5o  = exp  (ji+zF.o)  (4.5) 

where  F'  = 1 — (0.5)1/p.  For  p = 1 (the  lognormal  distribution),  zF>  = 0 and  t0  50  = exp (p). 

Density.  The  power-lognormal  probability  density  is  the  derivative  of  the  cdf  (4.2);  namely, 

g(f,p,a,p)  = (pita)  «{[ln(0-M]/o}  W-tlnW-rf/ojr1;  (46) 

here  4>(z)  = (2ir)“1/2  exp(-z2/2)  is  the  standard  normal  probability  density.  Figure  4.3  displays  such 

densities  as  a function  of  r = r/r,  where  r = exp  (p).  The  densities  have  a variety  of  shapes  depending 
on  the  p and  a values.  These  densities  are  unimodal.  For  p=  1,  the  density  is  lognormal. 

Mode.  The  mode  of  a power-lognormal  distribution  is  the  value  tm  where  the  probability  density 
peaks.  It  is  the  most  probable  value.  tm  is  the  solution  of  dg(t;p,o,p)ldt  = 0;  namely,  the  solution 
zm  = [\n(tm)—p\lo  of 

(zm+a)$(-zj  = -0-1).  (4.7) 

Here  zm  is  a function  of  p and  a.  For  p = 1 (the  lognormal  distribution),  zm  = —a  and  tm  = exp(p— a2). 

Mean.  The  mean  of  the  power-lognormal  distribution  is 
ET  = | “f  g(t;p,o,p)  dt 

= | “(p/a)  <*>{[ln (r)-/x]/a)  [c^ -[ln(r)-M]/a}]p_1  dt.  (4-B) 

= r|°°  pew <fr(u){$(-u)}p~x  du  = rE{o,p). 

Here  r = exp (p)  is  the  scale  parameter,  and  E{o,p)  is  a function  of  a and  p but  not  p = ln(r).  David 
(1981,  p 288)  references  tables  of  E{a,p)  for  integer  p and  a = 1.  The  integral  can  be  evaluated  with 
numerical  quadrature. 

Variance  and  standard  deviation.  The  variance  of  the  power-lognormal  distribution  is 
V(T)  = g(t;p,a,p)  dt 

= j ~(t-ET)2(plot)  <*>{[ln(r)-,*]/a}  [<&( -Dn(f)-^/ff}]'>-1  dt,  (4-9) 

= r2l“je’u-E(a,p)]2p  m [4<-h)]'’-1  du  - r2V(a,p). 

Here  r = exp (p)  is  the  scale  parameter,  and  V(o,p ) is  a function  of  o and  p but  not  p=ln(r).  David 
(1981,  p288)  references  tables  of  the  integral  for  integer  p and  a=  1.  The  integral  can  be  evaluated  with 
numerical  quadrature.  The  standard  deviation  of  the  power-lognormal  distribution  is 
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Fig.  4.3  Power-lognormal  probability  densities. 


sd(T)  - [V(T)]U1  = t ma,p)]'n.  (4I0) 

Forp=l,  the  lognormal  distribution,  V(a,l)  = [exp((T)  — 1]  exp((T),  V(T)  = r2  [exp(<r)-l]  x expCa2), 
and 


sd{T)  = r[exp(a2)-l]1/2exp(a2/2). 


(4.11) 


is 


Hazard  function.  The  hazard  function  (or  "failure  rate")  of  the  power-lognormal  distribution 


h(t\iL,o,p ) = g(t;p,o,p)/G(r,p.,o,p ) = p-h(t;p,a,  1); 


(4.12) 


here  h(t;ix,o,l)  = {\lot)4>{[\n(t/f)\/o}  / ${-[ln(r/T)]/(7}  is  the  lognormal  hazard  function  and  r = exp(/r) 
is  the  scale  parameter.  Rewritten 

h(t‘,p,a,p)  = p-/2(r;0,a,l)  (4-13) 

where  r = t/r.  Figure  4.4  shows  h(r;0,a,l)  for  selected  values  of  o.  These  are  a multiple  of  the 
lognormal  hazard  functions.  Thus,  they  have  the  same  behavior  as  lognormal  hazard  functions.  In 
particular,  they  are  all  zero  at  t= 0,  increase  to  a maximum,  and  then  decrease  back  toward  zero  with 
increasing  time.  For  most  o values  seen  in  practice,  the  maximum  occurs  in  the  far  upper  tail.  Thus, 
for  most  practical  purposes,  the  failure  rate  increases  as  the  population  ages. 


raio-i 


SlflM”  0.2.  0.4,  0.7.  1.0 


Fig.  4.4  Power-lognormal  hazard  functions. 
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5.  MAXIMUM  LIKELIHOOD  THEORY  FOR  FITTING  THE  POWER-(LOG)NORMAL 

MODEL  TO  DATA 

5.1  Introduction 

Purpose.  This  section  presents  maximum  likelihood  (ML)  theory  for  fitting  the  power- 
(log)normal  model  to  data.  As  explained  in  detail  by  Nelson  (1982,  1990),  the  ML  method  has  good 
properties  including 

• It  is  versatile;  that  is,  it  applies  to  most  models  and  types  of  data  including  censored  and  interval  data. 

• It  has  good  statistical  properties;  namely,  the  method  yields  good  estimators  and  confidence  limits. 

Overview.  The  contents  of  this  section  are: 

• Specimen  and  sample  likelihoods  (Section  5.2), 

• Maximum  likelihood  estimates  (Section  5.3), 

• Fisher  and  covariance  matrices  (Section  5.4),  and 

• Confidence  limits  (Section  5.5). 


5.2  Likelihoods 

Introduction.  This  section  presents  specimen  and  sample  (log)  likelihoods  for  the  power- 
(log)normal  model  and  observed,  left  censored,  right  censored,  and  interval  data.  Terminology  for  life 
data  is  used,  but  the  model  may  be  fitted  to  strength  and  other  types  of  data.  Needed  background  is  in 
Section  3 on  the  power-normal  distribution.  As  background,  the  following  paragraph  describes  types  of 
sample  data. 

Data.  The  types  of  censored  life  data  are  defined  here.  For  specimen  i with  size  Sit  failure  may 
occur  at  an  exactly  observed  failure  age  (life)  yf.  If  an  age  tt  is  from  a power-lognormal  distribution,  then 
one  uses  the  log  age  yf  = ln(rf)  throughout.  If  the  specimen  has  run  a time  y,  without  failure,  then  the 
failure  age  is  beyond  and  is  censored  on  the  right  at  age  y,-;  then  yf  is  called  the  survival  age  without 
failure.  If  the  failure  occurs  at  some  unknown  age  before  a known  first  inspection  age  yf,  then  the  failure 
age  is  censored  on  the  left  at  age  yt.  If  the  failure  occurs  at  some  unknown  age  between  two  adjoining 
inspection  ages  (y,-  < yf),  then  the  failure  age  is  interval  censored  and  in  the  interval  (y„y/).  For  a 
specimen,  each  such  type  of  data  has  a corresponding  likelihood  for  the  power-(log)normal  model,  as 
described  in  the  following  paragraphs. 

Observed.  The  likelihood  L,  for  specimen  i of  size  St  that  has  an  observed  failure  at  (log)  age 
y,-  is  the  probability  density  (3.7);  namely, 

l, = (kS/o)  Wy-iiVo]  {4<-0’i-M)/«];'‘Srl.  (51) 

This  and  other  likelihoods  are  functions  of  the  datay,  and  S,  and  the  model  parameters  \i,  a , and  k = 1 IS0 
where  S0  is  the  size  of  specimens  with  a (log)  normal  distribution.  The  log  likelihood  is 

2,.  - in(L(.)  = ln(*S/<7)  - ln(2,r)1/2  - \[(yrli)la)2  - (*S,-l)fl[0<rfi)/<j];  (5.2) 
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here 


H{z)  - -In m-z)] 


(5.3) 


is  the  cumulative  hazard  function  of  the  standard  normal  distribution. 


.KS: 


Right  censored.  The  likelihood  L,  for  specimen  i of  size  St  that  is  right  censored  at  (log)  age  yit 
its  current  unfailed  age,  is  the  survival  probability  (3.1);  namely, 

(5.4) 

(5.5) 


The  log  likelihood  is 


^ = W-(yrn)/o]} 

% = -KS'UUy-ri/o]. 


Left  censored.  The  likelihood  for  specimen  i of  size  S,  which  is  left  censored  at  (log)  age  yh  its 
inspection  age  after  its  unknown  failure  age,  is  the  failure  probability  (3.2);  namely. 


The  log  likelihood  is 


L,  - l - {fl-ty-jO/a]}'* 


% =■  ln[l  - 


(5.6) 

(5.7) 


Interval  censored.  The  likelihood  L,  for  specimen  i of  size  St  which  is  known  to  have  failed  in 
the  interval  between  inspections  at  ages  yt  < y[  is  the  probability  of  failure  in  that  interval;  namely, 

L,  - (5-8) 

The  log  likelihood  is  = ln(L,-). 

Sample  likelihood.  The  specimen  likelihoods  for  a sample  of  n specimens  yield  the  sample 
likelihood  below.  The  sample  log  likelihood  yields  ML  estimates  and  confidence  limits  for  model 
parameters  and  other  quantities,  as  described  in  the  following  sections.  The  specimen  lifetimes  are 
assumed  to  be  a random  sample  of  statistically  independent  lifetimes  from  a power-(log)normal 
distribution.  The  sample  likelihood  is  the  probability  of  the  n specimen  outcomes  (observed,  censored, 
or  interval  data).  For  independent  specimen  lifetimes,  the  sample  likelihood  L is  the  product  of  the 
specimen  likelihoods;  namely, 

L = Lx  x 1^  x ...  x Ln.  (5.9) 

Here  the  sample  likelihood  is  a function  of  the  data  values  y,-  (and  y[)  and  St,  the  type  of  data,  and  the 
model  parameters  \i,  a,  and  k = l/50.  The  sample  log  likelihood  2?  = ln(L)  is  the  sum  of  the  specimen 
log  likelihoods;  namely, 

S£  = S£j  + + - + 9Pn.  (5.10) 
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5.3  Maximum  Likelihood  Estimates 


Introduction.  This  section  defines  maximum  likelihood  (ML)  estimates  and  presents  related 

theory. 


Definition.  The  maximum  likelihood  estimates  for  the  model  parameters  \i,  a,  and  k are  the 

A A A 

values  /x,  o,  and  k that  maximize  the  sample  likelihood  (5.9)  (or  equivalently  that  maximize  (5.10))  over 
the  range  of  mathematically  allowed  values  of  those  parameters.  For  many  models  and  types  of  data,  the 
maximum  likelihood  estimates  are  unique.  Also,  the  ML  estimates  (/x,a,/c)  have  an  asymptotic  joint 
normal  distribution  with  mean  vector  (/x,a,x),  the  true  population  values,  and  covariance  matrix  (5.20) 
below.  Nelson  (1982,  1990)  references  authors  who  state  regularity  conditions  that  the  model  and  data 
must  satisfy  to  assure  this.  For  some  models  and  types  of  data,  there  may  be  local  maxima;  then  the 
parameter  values  at  the  global  maximum  are  the  ML  estimates.  In  many  computer  programs,  the  ML 
estimates  are  iteratively  calculated  by  direct  numerical  search  for  the  maximum  of  (5.10);  Section  6 
describes  such  calculations.  An  alternative  calculation  involves  the  likelihood  equations  which  follow. 

Likelihood  equations.  For  some  models  and  types  of  data,  the  ML  estimates  can  be  found  by 
the  calculus  method.  That  is,  the  parameter  values  for  which  the  partial  derivatives  of  the  sample  log 
likelihood  with  respect  to  each  parameter  all  equal  zero  are  the  ML  estimates.  For  the  power-(log)normal 
model  and  complete  data  (all  specimens  have  exactly  observed  failure  ages  yf),  the  partial  derivatives  set 
equal  to  zero  are: 

o = ase%  = (i/a){Lz,.  + E(^f-1)/I(zf)}, 

o = aae /da  = (i  /o){-n  + Ez?  + E(*sri)z/z(z,.)},  (5-H> 

0 = d&ldK  = (n/K)  - ZSfliZf). 

Here  the  standardized  deviate  z,-  s=  (y.—  ^)la  is  a function  of  /x  and  a,  and  each  sum  runs  over  all  n 
specimens.  These  likelihood  equations  cannot  be  solved  explicity  for  /x,  a,  and  k.  They  must  be  solved 
iteratively. 

Expectations.  Mathematical  expectations  of  various  quantities  appear  in  the  ML  theory  below. 
Such  expectations  are  defined  here  for  complete  data;  Nelson  (1982,  1990)  defines  them  for  other  forms 
of  data.  Let  u(yt)  denote  a given  function  of  the  random  failure  age  y{  of  specimen  i.  Then  the 
expectation  of  the  random  quantity  u(y()  is 

E{u(y,)}  = J °a_Ji(yl)ftyi\W,K)  dy{,  (5-12) 

here./(  ) is  the  power-normal  probability  density  (3.8),  and  the  parameters  fi , a,  and  k have  the  true 
population  values. 

Relationships.  A useful  general  property  of  the  first  partial  derivatives  of  the  specimen  and 
sample  log  likelihoods  is  that  their  expectations  are  zero.  That  is, 

0 = (l/a)[E  E{Zj}  + L(jcSrl)£{/i(z,.)}], 

0 = E{dWdo}  = (1  /a)[-n  + lE{tf}  + £ («Srl)E{z/r(z,)}],  (5-13) 

0 = £{a9£/a*}  = (nU)  - Y.SjE{H(Zi)}. 
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In  these  formulas,  the  y{  are  regarded  as  random  quantities  and  the  parameters  are  equal  to  their  true 
(unknown)  population  values.  These  formulas  aid  in  the  evaluation  of  various  expectations  later.  For 
example,  for  a single  specimen  i,  the  third  (5.13)  equation  yields 

ms,)  - E{H(z,)}  » f " dy,  = (I/O)  f " dz,;  (5.14) 

J—oo  J — oo 

Here^jO,!,^,)  is  the  standard  power-normal  distribution  with  p,-  = kS{.  Similarly,  the  first  (5.13) 
equation  yields 

£Wz,)}  - £{z,}/(«S,)  = (k5  -1)-1(1/o)  f " z,  Az„0 .US,)*,.  (515) 

The  last  integral  is  the  mean  of  the  standard  power-normal  distribution  with  p(-  = kS(.  These  results  are 
used  to  evaluate  integrals  below. 

5.4  Fisher  and  Covariance  Matrices 


Introduction.  This  section  presents  the  theoretical  Fisher  information  matrix  and  the  theoretical 
asymptotic  covariance  matrix  of  the  ML  estimates.  The  subsequent  estimates  of  these  matrices  are  used 
to  obtain  approximate  confidence  limits  in  Section  5.5. 

Second  derivatives.  ML  theory  employs  the  second  partial  derivatives  of  the  sample  log 
likelihood.  For  the  power-normal  model  and  complete  data,  they  are: 

32a/a,x2  = (-l/ff2){n.5;(/(Srl)/!fe1.)[/i(z1)-z1.]}, 
dtydodn  = ( - 1 /<T2)  { 2 E z,  + E (K  ^ - 1 )A  (z,)  [ 1 [ A(z()  -z;]  ] } , 

a2a/3><dfi  = (SJa)h(z:), 

322/3<t2  = (-l/ff2){-/!+3Ez^+£(KSj-l)zj/t(zI-)[2+zJ.[/!(zi)-zi]]}, 

322/3*3ff  = (1  lo)£S,z,Ji(z,), 

SHf/S*2  = -Vk2. 

In  these  equations,  the  derivative  of  the  standard  normal  hazard  function  satisfies  h'iz,)  = /i(z,)[/2(z,)-z,]. 

Fisher  matrix.  The  theoretical  Fisher  ( information ) matrix  is  the  matrix  of  expectations  of  the 
negative  second  partial  derivatives;  namely, 


F = 


E{-d2<y>/dfi2} 

E{-d2g/dodn} 

E{-d2&/dKdn} 


Ei-d^/dudo} 

E{-d2&/do2} 

E{-d2Wdicdo} 


£{-a29£/dpax} 

E{-d2<£JdodK} 

E{-d2(£]dK2} 


(5.17) 


This  matrix  is  symmetric  and  a function  of  the  parameters,  which  equal  their  true  (unknown)  population 
values  here.  The  ML  estimate  F of  this  matrix  is  obtained  by  replacing  the  parameters  with  their  ML 
estimates  in  (5. 17).  Some  integrals  in  (5. 17)  have  not  been  evaluated  analytically.  For  the  ML  estimate, 
the  integrals  can  be  evaluated  by  numerical  quadrature  when  the  true  unknown  parameter  values  are 
replaced  by  their  ML  estimates.  For  censored  data,  general  formulas  for  such  integrals  appear  in  Nelson 
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(1982,  p.  378)  and  require  a censoring  time  for  each  specimen,  even  those  that  failed.  Such  a planned 
censoring  time  may  not  be  known  for  field  data.  The  theoretical  matrix  is  useful  for  theoretically 
evaluating  proposed  test  plans  where  all  censoring  times  are  specified.  Another  estimate  of  F that  does 
not  require  all  censoring  times  and  is  easier  to  calculate  is  the  local  estimate 


F*  = 


-d29£/3/i2  -d2Q/dfido  -d2Q/dndic 
-d2Q/dodn  -d 2Q/dcr2  -d2<sZ/dodic 
-tf&ldicdfi  -d2&/dKdo  -d2&ldn2 


(5.18) 


This  is  the  matrix  of  negative  second  partial  derivatives  ((5.16)  for  complete  data)  of  the  sample  log 
likelihood  where  the  caret  A indicates  that  the  parameters  are  set  equal  to  their  ML  estimates.  For 
example, 

-d2&ldKdfi  = -(1  loyES/iiZ,)  (5-19) 

where  z{  = (yi—^/a.  Nelson  (1982)  gives  other  (equivalent)  formulas  for  the  theoretical  Fisher 
information  matrix  and  other  estimates  for  it.  This  local  estimate  of  the  Fisher  matrix  does  not  require 
a censoring  time  for  each  specimen,  only  for  the  censored  specimens.  Thus  it  can  always  be  used,  even 
when  the  ML  estimate  F cannot  due  to  unknown  censoring  times  for  failed  specimens. 


Covariance  matrix.  The  theoretical  asymptotic  covariance  matrix  2 of  the  ML  estimates  \l,  a, 
and  k is  the  inverse  of  the  theoretical  Fisher  information  matrix;  namely, 


z 


Var(/t)  Cov(ji,d ) Cov(/t,x) 
Cov(a,/t)  Var(a)  Co v(d,k) 
Cov(x,/i)  Co v(k,a)  Var(x) 


= F'1 


(5.20) 


Nelson  (1982,  pp.  384-386)  gives  a heuristic  argument  for  this  result.  2 is  a function  of  the  parameters 
H,  a,  and  k set  equal  to  their  true  (unknown)  population  values.  The  ML  estimate  of  the  covariance 
matrix  is  the  inverse  of  the  ML  estimate  (5.17)  of  the  Fisher  matrix;  namely,  2 = F-1.  Similarly,  the 
local  estimate  of  the  covariance  matrix  is  the  inverse  of  the  local  estimate  (5.18)  of  the  Fisher  matrix; 
namely,  2*  = (F*)-1.  An  estimate  of  the  covariance  matrix  yields  approximate  confidence  limits  for 
parameters  and  functions  of  parameters  as  described  in  Section  5.5. 


5.5  Confidence  Limits 

Introduction.  This  section  presents  approximate  confidence  limits  for  model  parameters  and 
functions  of  them. 

Function  estimate.  Often  one  wishes  to  estimate  some  function  6 = d(ji,a,K ) of  the  model 
parameters.  For  example,  the  F fractile  of  the  power-normal  model  for  specimens  of  size  S is  the 
function  tjf  = /t  + z ra  where  F = 1 — (1  — F)VkS  and  zF  = <i>-1(F').  Also,  for  example,  the  "normal 
size"  is  the  function  S0  = 1 Ik  of  just  k.  The  ML  estimate  of  the  true  function  value  d = 6(ji,o,k ) is  the 
function  evaluated  at  the  ML  estimates  of  the  parameters;  namely,  6 = 0(ll,o,k).  For  example, 
F'  = 1 — (1  -F)1  K and  tjf  = n + Under  certain  regularity  conditions  that  hold  for  many  models 
and  data,  the  asymptotic  cumulative  distribution  of  6 for  "large"  samples  is  close  to  a normal  cumulative 
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distribution  with  mean  6 (the  true  population  value)  and  variance  Var(0)  given  below  in  (5.22).  This 
asymptotic  normal  distribution  yields  the  approximate  confidence  limits  in  the  next  paragraph. 

Confidence  limits.  Based  on  the  asymptotic  normal  distribution  of  6,  lower  and  upper  two-sided 
approximate  100C%  normal  confidence  limits  for  the  true  value  6 are 

6 =0  - zc.[Var(#)]I/2,  6 = $ +zc-[V ar(0)]1/2;  (5.21) 

here  zc ■ is  the  standard  normal  C'  = (1 + Q/2  fractile,  and  Var(0)  appears  below  in  (5.22).  In  practice, 
one  uses  in  (5.21)  an  estimate  of  Var(0)  given  below. 

Variance.  This  paragraph  presents  the  formula  for  the  theoretical  asymptotic  variance  Var(0). 
It  is  assumed  that  the  function  6(ji,o,k ) has  continuous  partial  derivatives  with  respect  to  the  parameters 
fi,  a,  and  k.  Then 


Var($)  = [dd/dfi  36 Ida  3613k}  2 

= (36l3ii)2\2s{p.)  + 2 (3d/3iJ.)(dd/do)Co\(jL,d) 

+ (36  l3a)2Vai(o)  + 2 (36  ldo)(36  /3k)Co\(o  ,k) 

+ ( 36/3K)2Vai(k ) + 2 (36 I3k)(36 l3v)Co\(iL,k). 

Here  the  partial  derivatives  and  variances  and  covariances  are  all  evaluated  at  the  true  values  of  the 
parameters.  The  local  (or  maximum  likelihood)  estimate  of  this  variance  is  given  by  (5.22)  where  the 
partial  derivatives  are  evaluated  at  the  ML  estimates  of  the  parameters  and  the  variances  and  covariances 
are  replaced  by  their  local  (or  maximum  likelihood)  estimates.  In  practice,  the  local  estimate  is  usually 
easier  to  calculate.  Either  estimate  is  used  to  obtain  the  confidence  limits  (5.21). 

Covariance.  Suppose  that  one  wishes  to  estimate  another  function  \ p = i J/(jjl,o,k)  of  the  model 
parameters.  Its  ML  estimate  \j/  = \(/(Jl,o,k ) and  6 = 6(Jl,o,k)  have  an  asymptotic  joint  normal  cumulative 
distribution  for  large  samples.  That  joint  distribution  has  a mean  vector  (6,\J/),  the  true  population  values, 
variances  given  by  (5.22),  and  covariance 

Cov(M)  = (36/dii)(dt/dn)Var(ti)+(d6ldo)(d\J//da)Vai(d)+  (36  I3k)(3\PI3k)W2s(k) 

+ [(36 ldn)(3\J//do)+(dd  /da)(3\J//dn)]  Cov(/t,a)  (5  23) 

+ [(36 /3o)(d\l/ldK)+(36 ldK)(d\f//do)\  Co\(o,k) 

+ [(36 ldii)(d\p/dK)+(dd /dK)(d\p/dn)\  Cov(ji,k). 

This  result  is  used  in  Section  6 for  transformed  parameters  used  for  convenience  in  place  of  the  natural 
model  parameters  n,  a,  and  k. 

LR  limits.  The  approximate  normal  limits  (5.21)  are  simple  to  calculate  for  most  models  and 
types  of  data,  and  they  are  widely  used.  In  practice,  such  limits  usually  enclose  the  corresponding  true 
value  with  a probability  below  the  stated  confidence  100C%.  Calculation  of  likelihood  ratio  (LR) 
intervals  is  described  by  Nelson  (1990,  pp.  297-301).  Such  intervals  usually  have  a confidence  level 


36 1 dii 
36  Ido 
3613k 


(5.22) 
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closer  to  the  stated  100 C%.  However,  their  computation  is  more  complex,  and  they  are  not  yet  available 
in  POWNOR  and  in  most  computer  programs  for  ML  fitting. 

Improved  normal  limits.  The  approximate  normal  confidence  limits  (5.21)  can  often  be 
improved  by  working  with  a transformed  value  f = £(0)  and  treating  its  ML  estimate  f = £(0)  as 
normally  distributed.  Then  one  calculates  the  limits 

f = f - zc.[Var(r)]1'2.  r - f + zc.[Var(r)]1/2  (5.24) 

Here  Vard")  = (df/30)2  Var(0)  and  is  estimated  as  described  above.  Then  the  limits  for  0 are  obtained 
by  inverting  the  monotone  function  £(  ) to  get  0 =r~l(f ) and  d =f_1(r) , if  is  a*1  increasing  function 
of  0.  Otherwise,  0 =r"1(f)  and  d=£~l(£),  if~f(0)  is  a decreasing  function  of  0.  For  example,  in  the 
power-(log)normal 'model,  a useful  transformation  is  f = ln(or).  This  transformation  assures  that  the 
confidence  limits  for  a are  positive,  whereas,  (5.21)  may  yield  a negative  lower  limit,  which  is 
unacceptable  since  a > 0.  Such  transformed  limits  may  have  a confidence  level  that  is  closer  to  the 
stated  100C%. 

6.  NUMERICAL  METHODS 

Introduction.  This  section  describes  numerical  methods  that  facilitate  the  calculations  involved 
in  the  ML  theory  in  Section  5.  Many  formulas  and  ideas  in  Section  5 are  computationally  naive  and  must 
be  altered  as  shown  here.  In  particular,  finding  the  ML  estimates  by  iteratively  maximizing  the  log 
likelihood  is  difficult  unless  the  model  is  reparametrized.  Even  the  best  optimization  routine  has  difficulty 
optimizing  the  log  likelihood  of  a poorly  parametrized  model,  and  the  natural  model  parameters  are  often 
a poor  choice.  As  background,  Ross  (1990)  discusses  such  numerical  problems  and  solutions  in  detail, 
and  Nelson  (1982,  pp.  386-395)  provides  further  solutions.  This  section  will  interest  only  a few  who  are 
concerned  with  numerical  methods. 

Starting  values.  The  following  method  uses  percentile  estimation  to  obtain  starting  values  n0, 
o0,  and  k0  for  the  ML  iteration.  The  method  employs  three  selected  sample  order  statistics  y,-  < yr  < 
yr  with  corresponding  probability  plotting  positions  P < P'  < P".  The  next  paragraph  deals  with  the 
choice  of  i,  i',  and  i" . For  a complete  sample  of  n specimens  of  size  S,  the  plotting  position  of  order 

statistic  i is  Pt  = ( i-0.5)/n . The  P fractile  of  the  power-normal  distribution  is 

tjp  = n + zra  (6.1) 

where  ir  = 1— <21/<s  and  Q — 1 —P.  For  specified 

Q = \-P  = (n-i+0.5)/n,  Q'  = 1-P'  = (n-i '+0.5)/n,  Q~  = 1-P'  = (n-i~+0.5)/n,  (6-2) 

the  differences  of  the  corresponding  fractiles  are 

VP~-  VP-  = (V  V)a’  Vp  - VP  = (V  zt)° • (6-3^ 

The  ratio  of  these  equations  yields 

(Vp-  - Tip-)  / (t Ip  - Tip)  = [Z^qJ/kS  - Z^qI/kS]  / [Z^qJ/kS  - Zx_^-ksY  (6-4) 

Substitute  the  order  statistics  y„  yr,  and  yr  for  the  corresponding  fractiles  7jP,  tjp.,  and  7 iP»  on  the  left 

side  of  this  equation,  and  iterate  on  b = 1 I(kS)  on  the  right  side  to  find  the  value  b0  that  satisfies  this 
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Then  use  k0  = 1 /(boS)  as  the  starting  value.  Then  substituting  the  order  statistics  j,  and  yr  for  t]p  and 
rip.  into  (6.3)  yields  a starting  value  60;  namely, 

°o  =0v  ~ ~ zt>-  (6-5) 

where  ir  = l-^°and  i V = \-Q'b° . Finally,  substituting  the  order  statistic  yt  for  tjp  in  (6.1)  yields 

£0  = y>  " zfdo-  (6-6) 

These  starting  values  are  the  percentile  estimates  based  on  the  order  statistics  yh  yr,  and  yr.  These 
starting  values  can  be  used  to  calculate  starting  values  for  another  parametrization  of  the  model.  This 
method  extends  readily  to  multiply  (right)  censored  samples,  but  must  be  adapted  to  interval  data.  Treat 
each  failure  in  an  interval  as  if  it  occurred  at  the  midpoint  of  its  interval,  and  use  the  method  above.  This 
yields  crude  starting  values. 

Order  statistics.  For  a complete  sample,  the  order  statistics  used  in  the  previous  paragraph  are 
chosen  as  follows.  For  a (log)normal  distribution  (p=  1),  the  asymptotically  "best"  (as  n -*  00)  two  order 
statistics  for  estimating  a have  ranks  i and  i " satisfying 

0.0694  = (i-0.5)/n,  0.9306  = (*'-0.5)//2.  (6-7) 

The  solutions  are 

i = 0.0694/2+0.5,  r = 0.9306/2+0.5,  (6-8) 

which  each  must  be  rounded  to  the  nearest  integer.  Use  the  rounded  values  in  (6.8)  to  calculate  the 
plotting  positions  in  the  previous  paragraph.  For  the  middle  order  statistic  use  the  sample  median;  that 
is,  2'  = nil  for  even  n and  V = (/2+l)/2  for  odd  n,  and  P'  = Q'  = 0.50.  For  a (multiply)  censored 
sample,  calculate  plotting  positions  for  the  failure  times  as  described  by  Nelson  (1982,  Chap.  4),  and  use 
one  of  the  early  failures  Pt,  one  of  the  last  failures  Pr  and  one  about  midway  between  them,  Pv  — (P,- 
+ Pr)l  2. 


Iteration  parameters.  The  theory  in  Sections  3,  4,  and  5 employs  the  natural  model  parameters 
p,  a,  and  p (or  equivalently  S0  or  k = l/50).  These  parameters  are  mathematically  convenient  and  easy 
to  interpret,  especially  p and  a,  which  are  well  understood  as  parameters  of  the  normal  distribution 
(p=  1).  For  purposes  of  calculating  ML  estimates  by  optimizing  a sample  log  likelihood,  other  parameters 
work  better.  POWNOR  uses  the  parameters  p,  5 = ln(a),  and  u>  = ln(50)  in  the  log  likelihood.  That 
is,  exp(S)  replaces  a,  and  exp(oo)  replaces  S0.  The  new  parameters  all  range  from  - 00  to  + 00 . Thus 
the  optimization  with  respect  to  p,  5,  and  w is  unconstrained;  whereas  iterated  values  of  a and  S0  would 
have  to  be  constrained  to  be  positive,  which  complicates  the  optimization  calculations.  Moreover, 
confidence  limits  (5,6)  and  (w, w)  convert  into  positive  confidence  limits  a = exp(5),  <7  = exp(5), 
So  = exp(w),  and  So  = exp(ai)  avoiding  negative  lower  limits.  Also,  for  various  models,  such 
transformations  tend  to  give  actual  confidence  levels  that  are  closer  to  the  intended  levels  than 
untransformed  parameters  (Nelson  (1982,  p.393)). 

Optimization  of  the  log  likelihood.  The  ML  estimates  of  the  model  parameters  are  the  parameter 
values  that  maximize  the  sample  log  likelihood,  which  is  calculated  as  described  in  Section  5.2.  The 
program  employs  the  Powell  (1964)  method  to  iteratively  search  for  the  values  of  the  model  parameters 
that  maximize  the  log  likelihood.  This  method  does  not  use  the  explicit  derivatives  of  the  log  likelihood 
function  taken  with  respect  to  the  parameters.  The  user  must  provide  the  Powell  method  with  starting 
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values  for  the  ML  estimates  of  model  parameters.  A method  to  obtain  starting  values  for  the  ML  iteration 
is  given  above. 

Trimmed  standardized  deviates.  In  searching  for  the  parameters  that  maximize  the  log 
likelihood  function,  POWNOR  works  with  standardized  deviates  zt  = (y-folo  of  observations.  If  the 
starting  values  of  the  model  parameters  are  far  from  their  ML  estimates,  the  standardized  deviates  are 
large.  Such  large  values  can  cause  overflow  problems  and  slow  down  or  prevent  the  convergence  of  the 
iterative  method.  POWNOR  avoids  this  through  the  use  of  trimmed  standardized  deviates,  as  described 
in  Nelson  (1982,  p.  394).  POWNOR  uses  the  trimmed  standardized  deviates  until  the  iteration  is  close 
to  converging.  Then  the  program  switches  to  using  the  (untrimmed)  standardized  deviates  of  the 
observations  and  completes  the  ML  iteration. 

Standard  normal  cdf  and  its  inverse.  As  described  in  Section  5,  the  log  likelihood  function  for 
a censored  value  contains  the  standard  normal  cdf.  POWNOR  uses  the  double-precision  IMSL  (1987) 
subroutine  DNORDF  to  evaluate  the  standard  normal  cdf. 

Fisher  matrix.  After  POWNOR  obtains  the  ML  estimates,  it  calculates  the  local  estimate  of  the 
Fisher  (information)  matrix.  The  elements  of  this  matrix  are  the  negative  second  partial  derivatives  of 
the  log  likelihood  taken  with  respect  to  the  (transformed)  model  parameters  and  evaluated  at  their  ML 
estimates.  POWNOR  uses  the  STATLIB  subroutine  HESS  which  obtains  these  derivatives  numerically 
by  a perturbation  calculation  using  only  the  values  of  the  log  likelihood  function. 

Covariance  matrix.  The  local  estimate  of  the  covariance  matrix  for  the  ML  estimates  of  the 
model  parameters  is  the  inverse  of  the  local  estimate  of  the  Fisher  matrix.  The  program  uses  the  IMSL 
(1987)  subroutine  DLINDS  to  invert  the  local  estimate  of  the  Fisher  matrix. 

Confidence  limits.  The  program  calculates  approximate  confidence  limits  for  n,  ln(a),  and  ln(S0), 
using  their  ML  estimates  and  local  estimates  of  their  variances.  This  assures  positive  lower  limits  for  a 
and  S0.  Calculation  of  confidence  limits  for  a distribution  percentile  requires  the  inverse  of  the  standard 
normal  cdf.  The  program  uses  the  IMSL  (1987)  subroutine  DNORIN  to  evaluate  the  inverse  of  the 
standard  normal  cdf. 

Reparametrization.  Numerical  maximization  of  the  log  likelihood  function  is  typically  facilitated 
by  orthogonality  of  the  model  parameters.  Two  parameters  are  defined  to  be  orthogonal  if  the 
(asymptotic)  covariance  of  their  ML  estimators  is  zero  (Cox  and  Reid  (1987)).  Ross  (1990)  describes  how 
to  use  the  information  matrix  to  obtain  (approximately)  orthogonal  parameters.  This  reparametrization 
multiplies  the  original  parameter  vector  by  the  upper  triangular  Cholesky  decomposition  of  their 
information  matrix.  The  routine  POWNOR_ORTHO  uses  the  information  matrix  for  fi,  ln(u),  and  ln(S0) 
evaluated  at  the  current  values  of  these  parameters  in  such  reparametrization.  If  the  current  values  of  the 
parameters  are  already  close  to  the  ML  estimates,  the  transformed  parameters  are  then  approximately 
orthogonal.  The  routine  POWNOR_ORTHO  usually  maximizes  the  log  likelihood  more  quickly  with 
respect  to  the  transformed  parameters. 
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